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Observability of half-quantum vortices and skyrmions in p-wave superconductors is an outstand- 
ing open question. Under the most common conditions, fractional flux vortices are not thermody- 
namically stable in bulk samples. Here we show that in chiral p-wave superconductors, there is a 
regime where, in contrast lattices of integer flux vortices are not thermodynamically stable. Instead 
skyrmions made of spatially separated half-quantum vortices are the topological defects produced 
by an applied external field. 



PACS numbers: 74.20.Rp, 74.25. Dw, 74.25. Ha 

Higher broken symmetries in p-wave superconductors 
have inspired long-standing interest to realize topologi- 
cal defects more complicated than vortices. Much of the 
early discussions of various complex topological defects 
were in the context of superfluid 3 He. 1 Recently attention 
to these questions has raised dramatically in connection 
with superconductors which are argued to have p-wave 
pairing, such as Sr 2 Ru04 . The highly interesting possi- 
bility there, is connected with half-quantum vortices. 2 8 
Their statistics is non-Abelian and they could potentially 
be used for quantum computations. 9 Other kinds of topo- 
logical defects discussed in connection with spin-triplet 
superconductors are skyrmions 10 and hopfions. 11 In su- 
perconducting materials, the creation of these topological 
excitations is highly nontrivial. Superconducting com- 
ponents are coupled by a gauge field and there are also 
symmetry-reducing inter-component interactions. As a 
consequence fractional vortices have logarithmically or 
linearly divergent energies (see e.g. Ref. 8), while in- 
teger flux vortices have finite energy per unit length. 
Consequently, under usual conditions, half-quantum vor- 
tices are thermodynamically unstable in bulk systems. 
It was argued that complex setups, such as mesoscopic 
samples, are needed for their creation. 2,8,12 Recently it 
was claimed that a half-quantum vortex was observed 
in a mesoscopic sample of S^RuO^ 2 Other proposed 
routes to observe fractional vortices, invoke (i) thermal 
deconfinement, 3,6,13 (ii) potential materials with strongly 
reduced spin stiffness 4 and (iii) regimes very close to the 
upper critical magnetic field, where gauge-field mediated 
half-quantum vortex confinement is weak. 5 In some more 
general systems it was shown that fractional vortices 
could be thermodynamically stable near boundaries. 14 
Today the conditions under which half-quantum vortices 
and skyrmions 10 could be experimentally created in bulk 
superconductors still remains an outstanding open ques- 
tion. 

In this work we investigate the magnetic response of 
the Ginzburg-Landau model the has been widely ap- 
plied to Sr2Ru04. 15,16 Our considerations apply to two- 
dimensional systems or three-dimensional problems with 
translation invariance along the z-direction. Then the 
free energy density reads 



Figure 1. (Color on-line) - Numerically calculated texture of 
the pseudo-spin vector for a skyrmion carrying with a topo- 
logical charge Q = 2. As can be seen in the picture the 
skyrmionic topological charge density is confined in a closed 
domain wall. 

F(i> a ,A) = \V xA\ 2 (la) 
+ |A^>i| 2 + \D y ^ 2 \ 2 +7\D y ^\ 2 + 7 |A^2| 2 
+ 2 7 Re [(A^i)*A/^2 + (Dy^y D x ip 2 ] (lb) 

+ (2 7 -l)|^| 2 |^| 2 + ]T -i^ + ^x (lc) 

a=l,2 

+ Tl^i| 2 |^2| 2 cos(2(^ 2 - (^i)) . (Id) 

The different components of the order parameter are de- 
noted ^1,2 = |^i,2 \e i(f>1 > 2 ] D = V + ieA. The p-wave 
state is described here by a doublet of complex fields sub- 
jected to the the following symmetry breaking coupling 
: Re(^* 2 ^2) = |^i| 2 |^2| 2 cos(2((^ 2 - ipx)). The ground 
state breaks the U(l) x %2 symmetry, since the ground 
state phase difference is either 7r/2 or 37r/2. Gradient 
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Figure 2. (Color on-line) - A t her mo dynamically sta- 
ble skyrmion carrying two flux quanta, with e = 0.8 and 
7 = 0.5. Displayed quantities are, magnetic flux (A), the 
(inverted) energy density (B) and the sine of the phase dif- 
ference sin((^2 — <Pi) (C). On the second line, the densities 
of superconducting order parameter components \ipi\ 2 (D), 
\ip 2 \ 2 (E), and the 'doubled phase difference' Im(^! 2 ^1) (F). 
Panels (G) and (resp. H) on the third line are the supercur- 
rents associated with each component ipi (resp. ^2) of the 
order parameter (see Appendix A for definition). The last 
panel (I) shows the total supercurrent. 



terms (lb) make this model clearly anisotropic in the 
x2/-plane. The coefficient 7, controlling the anisotropy, 
should be 7 > 1/3 for S^RuC^ , according to. 15 The cou- 
pling constant e is a convenient quantity to parametrize 
the penetration depth of the magnetic field. The discrete 
%2 symmetry dictates that the system allows domain wall 
solutions interpolating between two regions with differ- 
ent phase-locking. Such domain walls are energetically 
expensive and thus not intrinsically stable. It was sug- 
gested that they could be observable if pinned by crys- 
talline defects. 17 Also domain walls formed as dynamic 
excitations inside vortex lattices were studies extensively 
in. 18 They could be experimentally observable in these 
setups since they pin half-quantum vortices. 17 ' 18 

Returning to the discussion of vortices one can ob- 
serve that the system (1) has U(l) x %2 broken sym- 
metry. Thus a single half-quantum vortex (with winding 
only in one of the phases) has linearly diverging energy 
and thus is not thermodynamically stable. 8 Also from 
this broken symmetry, the existence of skyrmionic ex- 
citations would not follow. The previous works required 
higher broken symmetry for the existence of skyrmions. 10 
However we show below that there is a considerable win- 
dow of parameters where the system (1) possesses what 



we term as a "skyrmionic phase". In that phase, mostly 
because of favorable competition between field gradients 
and potential and magnetic energies, the system does 
have thermodynamically stable skyrmions while ordinary 
integer flux vortex lattices are not thermodynamically 
stable. These skyrmions are bound states of spatially 
separated half-quantum vortices, connected by domain 
walls. Half-quantum vortices are linearly confined into 
integer vortices in a bulk sample because of the terms 
l^i | 2 1^2 1 2 cos(2((/?2 — <£i))- However on a (closed) domain 
wall, a composite vortex should split along this wall, since 
the above-mentioned term has there, unfavorable values 
of the phase difference. Indeed, such deconfining allows 
to reduce energetically unfavorable values of the phase 
differences. Because of this vortex splitting and resulting 
repulsive interactions, vortices trapped on domain wall 
can prevent the collapse of a closed domain wall. The 
main result of this paper is that we show that these ob- 
jects are characterized by an integer- valued skyrmionic 
topological charge and that they can be energetically 
cheaper than vortices. Such a skyrmion is displayed in 
Fig. 1, as a texture of a pseudo-spin vector field defined 
later on. 




Figure 3. (Color on-line) - A skyrmion carrying five flux 
quanta, with e = 0.8 and 7 = 0.4. Displayed quantities are 
the same as in Fig. 2, except panel (I) showing the gradient 
of the phase difference V<^i2, which is non zero at the domain 
wall. The skyrmion consists of ten spatially separated half- 
quantum vortices. It assumes a complicated non-symmetric 
structure due to a competition of a preferred geometry of a 
skyrmion with the anisotropics (lb). 

We investigated structures carrying N flux quanta 
(i.e. with each phase winding § Vip a = 2ttN) as functions 
of the gauge coupling e and the anisotropy parameter 7. 
Ground states, carrying a given number of magnetic flux 
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quanta, are computed numerically by minimizing the en- 
ergy within a finite element framework provided by the 
Freefem++ library. 19 See technical details in Appendix 
B. 

We find that if the penetration length is sufficiently 
large (i.e. at small values of the coupling constant e), 
the system indeed forms ordinary Abrikosov vortices in 
external field. On the other hand for sufficiently large e 
the system behaves as a type-I superconductor. However 
there is a regime in a wide range of intermediate cou- 
pling constants e, where integer flux vortices are more 
expensive than bound states of spatially separated half- 
quantum vortices connected by closed domain wall. Such 
configurations carrying different number of flux quanta 
are given in Figures 2, 3 and 4. The clearly visible 
preferred directions for supercurrents originate in the 
anisotropics (lb). The cores in different components do 
not coincide in space. This means fractionalization of 
vortices in this state. Each of the split cores carries a 
half of a flux quantum (for detailed calculations of frac- 
tional vortices flux quantization, see e.g. Ref. 8). 




Figure 4. (Color on-line) - A skyrmion with N = 8, e = 
0.6 and in the case of higher anisotropy 7 = 0.6. Displayed 
quantities are the same as in Fig. 2. 

The configurations found here are actually skyrmions, 
although it may not be obvious from the Figures 2, 3 and 
4. To prove that the solutions are skyrmions the two- 
component model (1) is mapped to an anisotropic non- 
linear a- model. 20 In that mapping the superconducting 
condensates are projected on the Pauli matrices a allow- 
ing to define the pseudo-spin vector n: 

n = {n x ,n y ,n z ) = where tf T = ^) . 

(2) 



The target space being a sphere, together with the one- 
point compactification of the plane defines the map n : 
S 2 S 2 . Such maps are classified by the homotopy class 
7T2(S 2 ) G Zj, so there exists an integer valued topological 
charge 

Q(n) = -— / n • d x ii x d y n dxdy . (3) 

For a skyrmion, Q = N, while Q = for ordinary vor- 
tices. The terms in (lc) and (Id), break the 0(3) sym- 
metry of the pseudo-spin n down to Z2. In a non-linear 
a-model, such anisotropy would undermine stability of 
the skyrmions. However this collapse does not occur in 
the Ginzburg-Landau model, because of the behaviour of 
the gradient energy, which is demonstrated below. 

The numerically computed topological charge (3) is 
found to be an integer (with a negligible relative er- 
ror of the order 10 -5 , due to the discretization) for the 
closed domain wall/vortex systems which are therefore 
skyrmions. The solutions shown in Figures 2, 3 and 4 
have skyrmionic topological charge Q = 2, Q = 5, Q = 8 
correspondingly. The terminology skyrmion is more in- 
tuitively obvious when the solutions are represented in 
terms of the pseudo-spin vector field n, as in Fig. 1. How- 
ever unlike skyrmions in non-linear a-model, here the 
skyrmionic topological charge density is mostly concen- 
trated on the half-quantum vortices and on the domain 
wall. 

The main result of this work is that skyrmions of the 
above type (and thus half-quantum vortices) can be less 
energetic than integer-flux ordinary vortices and ther- 
modynamically stable, in the chiral p-wave superconduc- 
tors. The critical external magnetic field H c \ for the 
formation of a flux-carrying topological defect is deter- 
mined by the condition where the Gibbs free energy 
G = Ed — 2 j B • H e dxdy becomes negative. Here Ed 
and B are the energy and magnetic field of the defect. 
H e denotes the applied field. Thus H c \ = Ed/2& where 

is the magnetic flux produced by the defect. The de- 
fects are thermodynamically stable if the critical external 
magnetic field's energy density H 2 X is smaller than the 
condensation energy. We investigated the energy depen- 
dence of the skyrmions on the number of enclosed flux 
quanta N. The energy of an integer flux vortex is used 
as a reference point. As shown in Fig. 5 panels (a) and 
(b), for low AT, the energy depends non-monotonically 
on N . This is because the preferred symmetry of small 
N configurations in some cases is in strong conflict with 
the anisotropics of the model. In the large- A" limit the 
energy per flux quantum gradually tends to some value. 
The main point here is that the energy per flux quan- 
tum for skyrmions is in certain cases smaller than that of 
vortices. This signals instability of vortex lattices with 
respect to skyrmion formation. 

Next, the thermodynamical stability of skyrmions is 
investigated. Results for N = 5 quanta are reported as 
a characteristic example, in Fig. 5 (c). We find that 
there are three regimes on the resulting phase diagram. 
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Figure 5. (Color on-line) - Upper panels show the depen- 
dence of the energy per flux quantum for skyrmions of dif- 
ferent topological charges Q (values are given in the units of 
the energy of one integer flux vortex). The N = 1 point at 
the origin corresponds to an ordinary vortex solution. Panel 
(a) shows calculations corresponding to different 7 for fixed 
e = 0.6, while (b) displays how the energy per flux quantum 
changes with e and N for fixed anisotropy parameter 7 = 0.7. 
The Q = 2 skyrmions are usually less energetically expensive 
than the Q = 3. This is because the Q = 2 skyrmions can 
be better aligned with the underlying anisotropics, than the 
Q = 3 skyrmions. 

The lower panel displays the phase diagram, calculated us- 
ing energy characteristics of Q — 5 skyrmions. The different 
colors refer to different physical properties. The type-I re- 
gion is shown by yellow shade. The lower part of the phase 
diagram shows regions where skyrmions (red) or vortex lat- 
tices (blue) form in applied external field. The phase diagram 
retains similar structure in calculations with different topolog- 
ical charges. With the increasing of the skyrmionic charge Q, 
the region where skyrmions are energetically preferred over 
vortex lattices grows slightly. These results apply either for 
two-dimensional systems or three dimensional systems with 
translational invariance along z-axis. In the latter case the 
energy should be understood as the energy per unit length of 
a skyrmion line (i.e. a skyrmion texture in xy plane which is 
invariant under translation along 2-axis). The discretization 
errors can be estimated by computing the total magnetic flux 
and comparing it to the exact value which follows from the 
quantization condition 2nN/e. This gives the relative accu- 
racy on the flux to be around 10 -5 . From that, the accuracy 
on the energy is estimated to be at least three order of mag- 
nitude smaller than the energy difference between skyrmions 
and vortices. 



When the penetration length is large {i.e. low e), the 
system shows usual the type-II superconductivity. When 
the penetration length is small, the system is a type-I 
superconductor. For intermediate values of the penetra- 
tion length, depending on the underlying anisotropies 7, 



Table I. Different contributions to the skyrmion energy per 
flux quantum. Q = 5 skyrmions are considered in this exam- 
ple. The results are compared with the contributions to the 
energy of a single vortex (which determines the lower bound 
on vortex lattice energy near the first critical magnetic field 
H c i). The gradient contribution E gra d is given by the in- 
tegrated (lb), the magnetic energy E mag by (la). The po- 
tential energy E pot is (lc) and E% 2 is (Id). First block, for 
which 7 = 0.8 and e = 0.4, corresponds to the state where 
skyrmions are thermodynamically stable but vortex lattices 
are not. Second block is for 7 = 0.6 and e = 0.2. It corre- 
sponds to a regime with standard Abrikosov vortex lattice. 
Here the skyrmions are local minima of the free energy func- 
tional. They are more expensive than vortices but, if formed, 
they are protected against decay by a finite energy barrier. 
In the second example the win in the kinetic energy is too 
small to overcome extra energy cost associated with domain 
wall formation and magnetic energy. 



the external field produces skyrmions rather than vor- 
tex lattices. To understand the instability of vortex lat- 
tices with respect to skyrmion formation, different con- 
tributions to energy are investigated in Table I. In the 
skyrmionic state, vortex lattice decay into skyrmions is 
driven by a win in gradient and potential energies al- 
though there is a loss in magnetic energy as well as the 
extra cost of producing a domain wall. 

The skyrmions we find are are structurally dif- 
ferent from skyrmions discussed in other kinds of 
superconductors 10 because of the different symmetry of 
the model. Another principal difference is the nature of 
the skyrmionic state, namely Ref. 10 proposed models 
where there are only skyrmionic solutions carrying two 
flux quanta. The latter forming stable lattices. In con- 
trast, the model we consider supports skyrmions with 
any integer value of topological charge. Importantly, the 
energy per flux quantum here is a sublinear function of 
the topological charge, which prohibits a ground state in 
the form of a lattice of the simplest skyrmions envisaged 
in Ref. 10. Instead our model predicts more complicated 
high-topological-charge skyrmionic structures. Also in 
type-II regime our model predicts metastable states of 
coexisting vortices and skyrmions. 

In conclusion we have shown that the phase diagram of 
chiral p-wave superconductors has a thermodynamically 
stable skyrmionic phase between type-I and the usual 
type-II regimes. This is despite the fact that the model 
has U(l) x Zj2 broken symmetry where naive symme- 
try arguments would rule out skyrmionic excitations. In 
the skyrmionic phase, the long sought-after half-quantum 
vortices acquire thermodynamic stability. These objects 
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can be detected with surface probes through their char- 
acteristic profile of magnetic field. The phase transi- 
tion into a skyrmionic state should be first order, be- 
cause the energy per flux quantum is decreasing with the 
skyrmionic topological charge. 

The possibly chiral superconductor S^RuC^ which is 
frequently described by the model (1) may have a pen- 
etration length which is slightly too large to fall into 
the skyrmionic phase. However in this case, the model 
predicts metastable skyrmionic excitations (which are 
slightly more energetic than vortices). Recently sporadic 
formation of objects with multiple flux quanta were re- 
ported in Fig. 2 of Ref. 21. Higher resolution scans of the 
magnetic field profile could confirm or rule out if the ob- 
served objects are skyrmions. Another scenario for flux 
clustering in this material is type- 1.5 superconductivity 



which can arise if to take into account its multi-band 
nature. 22 
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Appendix A: Technical details 

Functional variation of the free energy (1) with respect 
to the fields provides Euler-Lagrange equations of mo- 
tion. In particular, variation with respect to the gauge 
field defines the total supercurrent 

J = J(D +J (2). (A .i) 
The contribution of each component there, is given by 



4 1] = 


|lm [^(AcV'i ■+ 


-7^2)] 




4 1] = 




+ D x fo)] 




4 2) = 








7(2) _ 

Jy — 




■7^1)] • 


(A.2) 



The contributions \J^\ and |«7^| are displayed for the 
different solutions. Each contribution to the supercur- 
rent reflects the underlying anisotropics due to the com- 
plicated gradient terms (lb). 



Simple quantum vortex solutions 

In the main part of the text, skyrmionic excitations 
are discussed. For a comparison, we display in Fig. 6 
the solution for a usual integer quantum vortex solution 
in the model (1). Because of the anisotropics (lb) of 
the theory, the single quantum vortex is also non-axially 
symmetric. In contrast to vortices, skyrmions have a 
non-zero skyrmionic topological charge. Moreover, visual 
inspection of phase difference cf2 — (fi of the skyrmions 
provides further arguments of the qualitative difference 
from vortices. Indeed, for skyrmions, the phase difference 
(f2 — <Pi covers the range [— 7r/2, tt/2]. Thus it links both 
inequivalent ground states. Phase difference for integer 
vortex ranges only [0, tt/2]. 



ip± parametrization of the condensates 

Another parametrization of the condensate is some- 
times used in literature. Namely combinations 



V2 



(A.3) 



are used instead of The free energy functional (1) 

can be rewritten within this different representation (a 
special care is there needed in the redefinition of the pa- 
rameters). Let the total current be 

J = J (+) + J ( " } . (A.4) 
The contribution of each component in terms of the 




Figure 6. (Color on-line) - A single vortex solution, when 
e = 0.2 and the parameter 7 = 0.6. Displayed quantities are, 
magnetic flux (A), the (inverted) energy density (B) and the 
sine of the phase difference sin(<^2 — <£i) (C). On the second 
line, the densities of superconducting order parameter compo- 
nents |^i | 2 (D), |^2 1 2 (E), and the 'doubled phase difference' 
Im(V>* 2 V^) (F). Panels (G) and (resp. H) on the third line, 
are contribution \J^\ (resp. |«7^|) of ipi (resp. ^2) compo- 
nent to supercurrent (A.2). The last panel (I) shows the total 
supercurrent (A.l). The densities (D and E) and the super- 
currents (G and H) of the different components of an integer 
quantum vortex are very anisotropic. The anisotropics are 
less perceptible when considering the magnetic field (A) or 
the energy (B). The main difference with the skyrmions can 
be seen from panel C. Indeed the phase difference goes from 
zero at the vortex core to tv/2 faraway, instead of — tv/2 to 
7r/2 for a skyrmions. 



parametrization, is given by 



J: 



(1) 



J: 



(2) 



± -Re [^7(A^2 • 



J (±) 

Jy 



r(2) 



±-Re M(D y fo 



(A.5) 



where and are defined in (A.2). 

In order to compare the features of this new 
parametrization, Fig. 7 show the very same integer flux 
vortex solution as in Fig. 7, but using the parametriza- 
tion (A.3). In this new parametrization, the component 
-0+ (panel D) has zero ground state density, while ip- 
(panel E) recovers to non-zero ground state density far 
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Figure 7. (Color on-line) - A single vortex solution, for the 
same parameters as in Fig. 6. Quantities displayed now are 
using the ip± parametrization. Panel (A) is the magnetic flux, 
and (B) is the (inverted) energy density. The sine of the phase 
difference sin((^ + —<f-) (C). On the second line, the densities 
of superconducting order parameter components |^+| 2 (D), 
|^_| 2 (E), and the 'doubled phase difference' Im(^+ 2 ifr-) (F). 
Panels (G) and (resp. H) on the third line, are contribution 
| (resp. \J^\) of V>+ (resp. ip-) component to supercur- 

rent (A. 5). The last panel (I) shows the total supercurrent 
(A.4). 



from the vortex. Both ip+ and ip- components have a 
coinciding core singularity. Panel C, showing the sine 
of the phase difference sin(^_ — ^ + ) exhibit a non-zero 
winding around the vortex core. This kind of features re- 
produce the kind of single vortex solution obtained in. 16 
The skyrmion solutions are displayed within this field 
parametrization in Figures 8 and 9. They correspond to 
solutions provided in the main part. 



Appendix B: Numerical Methods — Finite element 
energy minimization 

We provide here details about the numerical methods 
which are used to construct skyrmions or vortex solu- 
tions of the two-components Ginzburg-Landau model (1). 
The variational problem is defined for numerical compu- 
tation using a finite element formulation provided by the 
Freefem++ library. 19 Discretization within finite element 
formulation is done via a (homogeneous) triangulation 
over Q, based on Delaunay-Voronoi algorithm. The do- 
main Q is chosen here to be a disc whose radius is much 
larger than the vortex/skyrmion size. In most cases the 
radius of the disc is 10 to 20 times larger than the size 




Figure 8. (Color on-line) - The thermodynamically stable 
skyrmion carrying five flux quanta, presented in Fig. 3 of the 
manuscript, using the ip± parametrization. Displayed here 
quantities are the same as in Fig. 7. 




Figure 9. (Color on-line) - The skyrmion carrying eight flux 
quanta, presented in Fig. 4 of the manuscript, using the ip± 
parametrization. Displayed here quantities are the same as in 
Fig. 7. 



of a single vortex. This guarantees that all the fields 
reach their ground state values at the boundary. This 
ensures that the topological solitons are not affected by 
the boundary. We performed the additional check that 
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solutions are not boundary artifacts by computing the 
energy on the boundary. When the algorithm converges, 
this quantity is of the order of the numerical accuracy, 
which indicates that the solutions do not interact with 
the boundaries of the numerical domain. Functions de- 
scribing the fields are decomposed on a continuous piece- 
wise quadratic basis over each triangle. The accuracy of 
such method is controlled through the number of trian- 
gles, (we typically used 3 ~ 6 x 10 4 ), the order of ex- 
pansion of the basis on each triangle (P2 elements being 



2nd order polynomial basis on each triangle), and also 
the order of the quadrature formula for the integral on 
the triangles. 

Once the problem is mathematically well defined, a nu- 
merical optimization algorithm is used to solve the varia- 
tional non-linear problem (i.e. to find the minima of J 7 ). 
Here we used a non-linear conjugate gradient method. 
The algorithm is iterated until relative variation of the 
norm of the gradient of the functional T with respect to 
all degrees of freedom is less than 10 -6 . 



